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1. INTRODUCTION 


I 

I 

In this paper, we describe and analyze numerical techniques that are | 

designed to approximate weak solutions of hyperbolic systems of conservation 
laws in several space dimensions. For sake of exposition, we shall describe 
these methods as they apply to the pure initial value problems (IVP) for a 
one-dimensional scalar conservation law 

(1.1) u^ + f(u)^ = 0, u(x,0) = 

To further simplify our presentation, we assume that the flux f(u) is a 

convex function, i.e., f"(u) > 0 and that the initial data uq(x) are i 

i 

piecewise smooth functions which are either periodic or of compact support. | 

Under these assumptions, no matter how smooth un is, the solution u(x,t) of | 

I 

the IVP (1.1) becomes discontinuous at some finite time t » t^,. In order to I 

extend the solution for t > t^,, we introduce the notion of weak solutions, 
which satisfy 

(1.2a) u dx + f(u(b,t)) - f(u(a,t)) - 0 

a 

for all b ^ a and t ^ 0. Relation (1.2a) implies that u(x,t) satisfies 
the PDE in (1.1) wherever it is smooth, and the Rankine-Hugoniot jump relation 

(1.2b) f(u(y + 0,t)) - f(u(y - 0,t)) = [u(y 0,t) - u(y - 0,t)] 

across curves x = y(t) of discontinuity. 


- 2 - 


It is well known that weak solutions are not uniquely determined by their 
initial data. To overcome this difficulty, we consider the IVP (1.1) to be 
the vanishing viscosity limit e + 0 of the parabolic problem 

(1.3a) (u®)^ + f(u®)^ = u^(x.O) = Uq(x), 

and identify the unique "physically relevant" weak solution of (1.1) by 

( 1.3b) u = lim , ^u® . 

e + u 

The limit solution (1.3) can be characterized by an inequality that the 
values u^ = u(y - 0,t), ® u(y + 0,t) and s = dy/dt have to satisfy; 

this inequality is called an entropy condition; admissible discontinuities are 
called shocks. When f(u) is convex, this inequality is equivalent to Lax's 
shock condition 


(1.4) 


a(u^) > s > a(Uj^) 


where a(u) = f'(u) is the characteristic speed (see [20] for more details). 
We turn now to describe finite difference approximations for the 


numerical solution of the IVP (1.1). Let Vj” denote the numerical 
approximation to u(xj,tjj) where xj = jh, tjj = nr; let vj^(x,t) be a 

globally defined numerical approximation associated with the discrete values 

{vj”},“<j<«, n>^0. 

The classical approach to the design of numerical methods for partial 


differential equations is to obtain a solvable set of equations for 
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by replacing derivatives in the PDE by appropriate discrete approximations. 
Therefore, there is a conceptual diffficulty in applying classical methods to 
compute solutions which may become discontinuous. Lax and Wendroff [21] 
overcame this difficulty by considering numerical approximations to 
the iveak ^OA/naCcLtion (1.2a) rather than to the PDE (1.1). For this purpose, 
they have introduced the notion of schemes in conservation form: 


(1.5a) 


■ ’j" - = ‘V ’“'j! 


here X = x/h and f^^ denotes 


(1.5b) 


- « f / 

'1+1/2 '''^1-k+r 


n \ 

... 


f (w^ , . . . ,w„, ) is a numerical flux function which is consistent with the 

X ^ IV 

flux f(u), in the sense that 


(1.5c) 


f (u,u, . . .u) = f (u) ; 


E^ denotes the numerical solution operator. Lax and Wendroff proved that if 
the numerical approximation converges boundedly almost everywhere to some 
function u, then u is a weak solution of (1.1), i.e., it satisfies the weak 
formulation (1.2a). Consequently discontinuities in the limit solution 
automatically satisfy the Rankine-Hugoniot relation (1.2b). We refer to this 
methodology as shock-capturing (a phrase coined by H. Lomax). 

In the following, we list the numerical flux function of various 3-point 
schemes (k = 1 in (1.5b)): 
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(i) The Lax-Friedrichs scheme [19] 


(1»6) f(Wj,W2) = Y + f(w2> " 7 (^2 - Wj)] 

(ii) Godunov's scheme [5] 


(1.7a) 


f(WpW2> = f(V(0;Wj ,W2)); 


here V(x/t; W|^,W 2 ) denotes the self-similar solution of the IVP (1 
the initial data 

I Wj X < 0 

. 

W 2 X > 0 

(iii) The Cole-Murman scheme [26]: 


(l»8a) f(Wj,W2> = Y [f(Wj) + f(w2> - |a(w^,W2> | (v /2 - w^] 


where 


(1.8b) 


a(Wj,W2) 


f(w 2 ) - f(w^) 

W 2 - Wj 


a(Wj) 


if Wj * ^2 


if Wj = W 2 


(iv) The Lax-Wendroff scheme [21]: 


1 W + Wo 

7(Wj,W2) “I ■ ^a(-^- 2 — ^)[f(w 2 ) - f(w^)]}. 


1) with 


(1.9) 


2 
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(v) MacCormack's scheme [24]: 


( 1 . 10 ) 


{^(^ 2 ) + f(Wj - X[f(w2> - f(Wj)])}. 


Let E(t) denote the evolution operator of the exact solution of (!•!) 
and let denote the numerical solution operator defined by the RHS of 
(1.5a). We say that the numerical scheme is r-th order accurate (in a 
pointwise sense) if its local truncation error satisfies 


( 1 . 11 ) 


E(t) • u - E, • u = 0(h^^^) 
n 


for all sufficiently smooth u; here t = 0(h). If r > 0, we say that the 
scheme is consistent. 

The schemes of Lax-Friedrichs (1.6), Godunov (1.7) and Cole-Murman (1.8) 
are first order accurate; the schemes of Lax-Wendroff (1.9) and MacCormack are 
second order accurate. 

We remark that the Lax~Wendroff theorem states that if the scheme is 
convergent, then the the limit solution satisfies the weak formulation (1.2b); 
however, it need not be the entropy solution of the problem (see [11]). It is 
easy to see that the schemes of Cole-Murman (1.8), Lax-Wendroff (1.9) and 
MacCormack (1.10) admit a stationary "expansion shock" (i.e., f(u^) = ^('^r) 
with ^(^L^ ^ ^ steady solution. This problem can be easily 
rectified by adding sufficient numerical dissipation to the scheme (see [25] 
and [10]). 

The cardinal problem that is yet to be resolved is the question of 
convergence of the numerical approximation. 
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2. LINEAR STABILITY AND CMfPllTATION OF WEAK SOLUTIONS 

Let us consider the constant coefficient case f(u) = au, a = const, in 
(1.1), l.e., 

(2.1a) + au^ = 0, u(x,0) = “^(x), 

the solution to which is 


(2.1b) 


u(x,t) = Uq(x - at). 


In this case, all the schemes mentioned in the previous section, (1.6) - 
(1.10), take the form 


( 2 . 2 ) 


n+1 

V. 

1 


k 


» I 

I*-k 


C V? 
i J-I 





> 


where are constants independent of j (C^^ are polynomial functions 

of the CFL number v = Xa). We note that in the constant coefficient case 
Godunov's scheme is identical to that of Cole-Murraan; the MacCormack scheme is 
identical to that of Lax-Wendroff . Since the numerical solution operator 
Ejj of these schemes in the constant coefficient case becomes a linear 
operator, we shall refer to these schemes as essentially linear or just 
"linear" schemes. 

Next we briefly review the convergence theory of linear schemes; we refer 
the reader to [29] for a detailed analysis. 

We say that the numerical scheme is stable if 
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(2.3a) 1 (Ej^)”ll < C for 0 £ nx ^ T, x = 0(h). 

The constant coefficient scheme (2.2) is stable if and only if it satisfies 
von Neumann's condition: 

(2.3b) 1^ C e"^^^|<l for all 0<5<it. 

il=-k ^ 

It is easy to verify that all the 3-point schemes (1.6) - (1.10) satisfy 
condition (2.3b) under the Courant-Friedrichs-Lewy (CFL) restriction 

(2.4) |v| = |Xa| < 1. 

and thus are linearly stable. The notion of stability (2.3a) is related to 
convergence through Lax^s equivalence theorem, which states that a consistent 
linear scheme is convergent if and only if it is stable. 

The accumulation of error in a computation with a linearly stable scheme 
(2.2) is linear, in the sense that if the local truncation error (1.11) is 
0(h^"^^), then after performing N = T/x = 0(h ^) time-steps, the error is 
O(h^), i.e., 

(2.5) u(Xj,Nx) - Vj = 0(h’^). 

An immense body of work has been done to find out whether stability of 
the constant coefficient scheme with respect to all ’^frozen coefficients" 
associated with the problem, implies convergence in the variable coefficient 
case and in the nonlinear case. 
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In the variable coefficient case, where the numerical solution operator 
is linear and Lax's equivalence theorem holds, it comes out that the stability 
of the variable coefficient scheme depends strongly on the dlssipativity of 
the constant coefficient one, i.e., on the particular way it damps the high- 
frequency components in the Fourier representation of the numerical solution. 

In the nonlinear case, under assumptions of sufficient smoothness of the 
PDE, its solution and the functional definition of the numerical scheme, 
Strang proved that linear stability of the first variation of the scheme 
implies its convergence; we refer the reader to [29] for more details. 

Although there is no rigorous theory to support the supposition that 
linearly stable schemes should converge in the case of discontinuous solutions 
of nonlinear problems, we find in practice that this is true in many (although 
not all) instances; when such a scheme fails to converge, we refer to this 
case as "nonlinear instability". The occurrence of a nonlinear instability is 
usually associated with insufficient numerical dissipation which triggers 
exponential growth of the high-frequency components of the numerical solution. 

Next we present two shock-tube calculations by the scheme (1.5) with the 
numerical flux 


(2.6) KwpW^) = j {f(w2) + f(Wj - X[f(w2) - f(wp]) 




The shock-tube problem is modelled by a Riemann IVP for the one-dimensional 
Euler equations of compressible gas: 


^t + 



X < 0 


(2.7a) 


u(x,0) = 


X > 0 
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where 

(2.7b) u = (p, qp, E)^, f(u) = qu + (0, p, qp)^ 

with 

(2.7c) P = (y - 1)(E - j pq^). 

Here p, q, p, and E are the density, velocity, pressure, and total 
energy, respectively. In these calculations, y = 1*^ and 

(2.7d) = (0.445, 0.3111, 8.928), = (0.5, 0., 1.4275). 

The exact solution to this Rlemann problem consists of a shock propagating to 
the right followed by a contact discontinuity and a left-propagating 
rarefaction wave; it is shown by a continuous line in Figures 1 and 2. The 
numerical solution of (2.6) is shown in Figures 1 and 2 by circles. 

Figure 1 shows the results of the second-order accurate MacCormack 
scheme, i.e., 0=0 in (2.6). Observe the large spurious oscillations at 

the shock and at the contact discontinuity — this is a Gibbs-like 
phenomenon. Note that although the rarefaction wave is computed rather 
accurately, there are some spurious oscillations at its right endpoint due to 
the discontinuity in the first derivative there. 

Figure 2 shows the results of the first-order accurate scheme (2.6) 
with 6=1. Observe that now the numerical solution is oscillation-free. 
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However, both the shock and the contact discontinuity are now smeared much 
more than the corresponding ones in Figure 1* Note the excessive rounding of 
the corners at the endpoints of the rarefaction wave. 

It is important to understand that the Gibbs-phenomenon by itself is not 
an instability; this is self-evident when we consider the constant coefficient 
problem (2.1) with discontinuous initial data Uq, However, in compressible 
gas calculations, where both density and pressure are restricted to have 
nonnegative values, the Gibbs phenomenon may cause the numerical solution to 
get out of the physical domain. Attempting to replace negative values of 
density and pressure by positive ones makes the scheme nonconservative and may 
result in an exponential growth of the solution. 

The comparison between Figure 1 (0=0) and Figure 2 (0=1) 

shows that the Gibbs phenomenon in the second-order accurate scheme can be 
controlled by the addition of a numerical viscosity term. To do so without 
losing the second-order accuracy. Lax and Wendroff [21] suggested to take in 
(2.6) 0 = 0 (Wj,W 2 ) of the form 

(2.8) eCw^.w^) = xlaCw^) - a(wpl; 

here a =» f'(u) and x is a dimensionless constant; observe that 
0=0 in the constant coefficient case. 

Numerical experiments showed that as x increases the size of the 
spurious oscillations decreases, but at the cost of increased smearing of the 
discontinuity. Furthermore, when x is fixed, the size of the spurious 

oscillations increases with the strength of the shock. These observations 
Indicate that the numerical viscosity term (2.8) does not have an approriate 
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functional dependence on the parameters that control the Gibbs phenomenon. 
Consequently, the choice of a suitable value of x is problem dependent, 
and the practical use of the numerical scheme requires several preliminary 
runs to "tune parameters". 

Ideally, we would like to have high-order accurate schemes that are 
capable of propagating a shock wave without having any spurious 

oscillations. In the scalar case, this can be accomplished by designing 
schemes to be raonotoniclty preserving, i.e., to satisfy 


(2.9) 


V 


monotone 



V monotone. 


Godunov [5] has considered this avenue of design in the constant coefficient 
case (2.1) and showed that monotonicity preserving ZindjCUt schemes (2.2) are 
necessarily only first order accurate. For some time this result has been 
perceived as saying that high-order schemes are necessarily oscillatory. Only 
much later was it realized that Godunov's result applies only to linear 
schemes and that it is possible to design noyLtLne.ctA high order accurate 
schemes that are monotonicity preserving (see [1], [22], [6], [23], [7], [2], 
and [30]). Schemes of this type are the "modern shock-capturing schemes" 
referred to in the title of this paper. 

In the rest of this paper we concentrate on the design and analysis of 
such highly nonlinear schemes. Even in the constant coefficient case these 
schemes are nonlinear to the extent that does not justify the use of local 
linear stability. Therefore, we shall start our journey into the nonlinear 
world by introducing the notion of total variation stability, which is more 
suitable to handle this type of schemes. 
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3. TOTAL VARIATION STABILITY AND TVD SCHEMES 

Glimm [4] has considered the numerical solution by a random choice method 
of an IVP for a system of conservation laws with initial data of small total 
variation, and proved existence of weak solutions by showing convergence of 
subsequences. Following ideas used in Glimm's convergence proof, we can 
formulate the following theorem for convergence to weak solutions. 

Theorem 3*1: Let v^ be a numerical solution of a conservative scheme 

(1.5). 


(i) If 


(3.1) 


TV(v^(.,t)) < C . TV(uq) 


where TV( ) denotes the total variation in x and C is a constant 
independent of h for 0 ^ t ^ T, then any refinement sequence h ^ 0 
with T = 0(h) has a convergent subsequence h^ -»• 0 that converges in 
to a weak solution of (1.1). 

(ii) If Vj^ is consistent with an entropy Inequality which Implies 
uniqueness of the IVP (1.1), then the scheme is convergent (i.e., all 
subsequences have the same limit, which is the unique entropy solution of the 
IVP (1.1)). 


We remark that unlike convergence theorems of classical numerical 
analysis, in which one shows that the distance between the solution and its 
numerical approximation vanishes as h 0, the convergence argument in the 
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above theorem relies on a combination of compactness and uniqueness; its 
relation to an existence proof is quite obvious (see [8] for more details). 

Next we demonstrate the use of Theorem 3.1 to prove convergence of 
schemes in conservation form (1.5) which are monotone, i.e., are of the form 


(3.2) 


1 u/n n \ nv 

') ■ >k> ■ • ” >j’ 


where H is a monotone nondecreasing function of each of its arguments in the 
interval [a,b], a = min v^, b = max Vj. We note that the schemes of Godunov 
(1.7), Lax-Friedrichs (1.6) and the first order scheme (2.6) with 9=1, 
are all monotone. 

We start by observing that the operator in (3.2) is order preserving 


(3.3a) 


u > V 




E. • u > E, • V. 
h — h 


Since Ey, is also conservative. 


(3.3b) I (E • v) = I V 

j ^ j ^ 

it follows then from a Lemma of Crandall and Tartar (see [3]) that Ej^ is 
f ^-contractive, i.e., for all u and v in 

(3.3c) HE^ • u - E. • vll, < Hu - vll. 

h hi,— t. 


(here SuS ~ I 

^1 j ^ 


Taking u to be a translate of v, i*e«, 
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u . 
J 


j+1 


for all j 


we get from (3.3c) that 

(3.4a) TV(E. • v) < TV(v) 

n — 

where 

(3.4b) TV(w) = I |w^^^ ■ ''jl* 

It follows then that the numerical solution satisfies (3.1) with C = 1; thus 
we have established the convergence of subsequences. To show that all limit 
solutions are the same, we can use an argument of Barbara Keyfitz in the 
appendix to [11], which shows that (3.3c) implies that the scheme is 
consistent with Oleinik's entropy condition. This shows that monotone schemes 
satisfy the requirements of Theorem 3.1 and thus are convergent. 

Unfortunately, monotone schemes are necessarily only first order accurate 
(see [11]). However, once we give up the requirement (3.3a) that Ejj be an 
order preserving operator and consider the larger class of schemes that 
satisfy only (3.4), i.e., schemes that are total-variation-diminishing (TVD), 
it becomes possible to obtain second order accuracy. Observe that TVD schemes 
are necessarily monotonicity preserving (see [7]). 

The following theorem provides an almost complete characterization of TVD 
schemes (see [7], [8], and [18]). 



-15- 


Theorem 3.2: Let be a numerical solution operator of the form 


k-1 


(3.5a) vj*' - vj + 


where 


(3.5b) 


A n n 

*1+1/2'' ■ ''l+l 


V 


n 

i’ 


and (j ) denotes some functional of evaluated at j . Then is 

TVD if (and only if)^ the following relations hold: 


(3.6a) C_j(j - 1) > C_ 2 (J - 2) ... > C_^(j - k) > 0 

(3.6b) > ••• > “Vl^j + k - 1) > 0 

(3.6c) ^-1^^ - 1) < 1. 


We turn now to consider the important case of k == 1 in (3.5), i.e., 


i 

i (3.7) 


n+1 


'j * ‘=-l‘l>*jtl/2’' 


>* 3 - 1 / 2 '' 


we refer to (3.7) as an essentially 3-point scheme, because the coefficients 
Cq( j ) and C_j(j) may depend on more than just v”, To see 


^ Theorem 3.2 is not a complete characterization of TVD schemes, since the 
representation of a given nonlinear scheme in the form (3.5) is not unique. 
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the relation between the form (3.7) and the conservation form (1.5) let us 
consider the scheme 

(3.8a) - vj - - r,.i/2) 

with 

^i+l/2 “ I ^^i ^i+1 " *^i+l/2^i+l/2^ 

It is easy to see that (3.8) can be rewritten in the form (3.7) with 

(3.9a) C_^(j) = j ^®j+l/2 ‘^j+1/2^ 

(3.9b) ^0^^^ " T ^^j-1/2 ~ *^j-l/2^’ 

here ®i+i/2 “ which is defined by (1.8b). 

Applying Theorem 3.2 to the scheme (3.8), we get that it is TVD if 

(3.10a) ^l^j+1/2^ — ^**j+l/2 — 

We turn now to outline the modified flux approach for the construction of 
second order accurate TVD schemes (see [7]). To simplify our presentation we 
choose in (3»10a) 


*^j+l/2 ~ l®j+l/2l ’ 


(3.10b) 
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thls makes (3.8) identical to the Cole-Murman scheme (1.8). We observe that 
the TVD property of this scheme does not depend on the particular value of 
f(u), but only on the CFL-like condition 


(3.10c) 


^l®j+l/2l - 


note that this condition involves only the g^d voMieJ> f j . Consequently, if 

mod 

we apply this scheme to a modified flux f™ = f j + gj , i.e., 


(3.11a) 


'j ^ ■ ^^^j+1/2 ■ ^j-1/2^ 


(3.11b) ~ 2 ^j+1 ®j+l ~ l®j+l/2 ^j+l/2l^j+l/2'^ ^ 


where 


(3.11c) ''^j+1/2 ~ ^®j+l ~ ®j^^^j+l/2^ ’ 

we can conclude that this scheme is TVD provided that 


(3.12) ^l®j+l/2 ^ ^j+l/2l - 

It is easy to verify by truncation error analysis that if 

g, » ha(a)u + O(h^) 

J * 


(3.13a) 
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where 

(3.13b) a(x) = Y lx| (1 - X|x|), 


then 

(3.13c) 


^j+1/2 “ ^j+1/2 


— LW 

where f is the numerical flux (1*9) of the second-order accurate Lax- 
Wendroff scheme. 

In [7] we have taken gj to be 


(3.14a) 






where 

(3.14b) m(x,y) 



min( |x| , |y| ) 
0 


if sgn(x) = sgn(y) = s 
otherwise 


Clearly g^ (3.14a) satisfies (3.13a) and consequently the resulting 

2 

scheme is second-order accurate, except at local extrema where the 0(h ) 
term in (3.13a) and (3.13c) fails to be Lipschitz continuous. 

Next we show that due to this particular definition of g^ , the modified 
flux scheme (3.11) which is second-order accurate, is also TVD under the 
original CFL restriction (3.10c); this follows immediately from the following 
lemma. 

Lemma 3.3. 


(3.15a) ( i) 


‘ <VV2’I 
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(3.15b) (il) X < 1 ^ X Yj+l/J < 1- 

Proof ; We note that m(x,y) (3.14b) satisfies |m(x,y)| < min ( |x|, |y|)« 


Consequently 


|gj+,-8jl < ISjl + l8j.ll < "i" ' IV'/zWi'’"!- IW2W2''"I> 

■f .dn ( 3/2.” I ) 


< 2 IVl/z” 


til 


which proves (3.15a). 

It follows therefore from (3.13b) and (3.15a) that 

X |a +y| < 1 l®l < lla| + 2X |a(a)| 

= x|a| + x|a| (1 -x|a|)< X|a| + 1 -xlaj = 1, 
which proves this lemma. 

We remark that the modified flux scheme (3.11), as the Cole-Murman scheme 
It Is derived from, admits a stationary "expansion shock" as a steady 
solution. Replacing 9j+i/2 “ l^j+l/Z^ (3.10b) by 


(3.16) 


*^j+l/2 “ “®’^^l®j+l/2l’ 


e > 0 
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results in a modified flux scheme which is entropy consistent (see [28]) and 
thus can be shown to be convergent by Theorem 3.1. 

The choice (3.14) of is by no means unique. It is easy to check 

that changing gj to be 


(3.14a)' 

with 

(3.14b)' 

or 

(3.14a)' 


8j - S VV^ > 

( X if |x| < |y| 

If |x| > |y| 

(> - X Aj^lyv") 


with 


(3.14b)" 


m(x,y,z) » m(m(x,y),z) , 


does not alter the relations (3.13a) and (3.15a) which makes the modified flux 
scheme (3.11) a second-order accurate TVD scheme, under the original CFL 
restriction (3.10c). 

The modified flux approach is not the only methodology to construct 
second order accurate TVD schemes (there are many ways to skin a nonlinear 
cat). In the next section, we shall describe the MUSCL scheme of van Leer 
[23]; other techniques are described in [30], [27], and [31]. Unfortunately, 
all TVD schemes, independent of their derivation, are only first order 
accurate at local extrema of the solution. Consequently, TVD schemes can be 
second-order accurate in the Lj sense, but only first order accurate in the 
maximum norm (see [14] for more details). 



- 21 - 


4. GODONOV-TYPE SCHEMES 

In this section we describe Godunov-type schemes which are an abstraction 
of Godunov's scheme (1.7) (see [5]) due to ideas in [23], [12], and [13]. 

We start with some notations: Let {1^} be a partition of the real line; 
let A(I) denote the interval-averaging (or "cell-averaging") operator 


(4.1) 


A(I) . w 


1 


/ w(y)dy ; 
I 


let Wj * denote w = {w^} . We denote the approximate 

reconstruction of w(x) from its given cell-averages {w^} by R ( x ;w). To 
be precise, R(x ; w) is a piecewlse-polynomial function of degree (r-1), 
which satisfies 


(4.2a) (1) R (x ; w) = w (x) + O(h^) wherever w is smooth 


(4.2b) (ii) A (Ij) • R ( • ; w) = w^ (conservation). 

Finally, we define Godunov-type schemes by 
(4.3a) = Adj”"^^) . E(x) • R(* ;v”) = (E^^ • v")^ 

(4.3b) v^° = % * 

here partition of the real line at time t^, and E(t) is the 

evolution operator of (1.1). 
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In the scalar case, both the cell-averaging operator 
solution operator E(t) are order-preserving, and consequently also total- 
variation diminishing (TVD); hence 

(4.4) TV (E. . w) < TV (R (• ; w)). 

n 

This shows that the total variation of the numerical solution of Godunov- 
type schemes is dominated by that of the reconstruction step. 

The original first-order accurate scheme of Godunov is (4.3) with the 
piecewise-constant reconstruction 

(4.5) R (x ; w) = Wj , for x c 1^ . 

Since the piecewise-constant reconstruction (4.5) is an order-preserving 
operation, it follows that E^ is likewise order preserving as a composition 
of 3 such operations; consequently the scheme is monotone. 

The second-order accurate MUSCL scheme of van Leer [23] is (4.3) with the 
piecewise-llnear reconstruction 


(4.6a) 


R(x; w) = w. + (x - y . ) • s. for x e I. , 

J 11 1 


where s. is defined by 
1 




(4.6b) 
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here y. denotes the center of I. . It is easy to verify that the particular 
•3 3 

form of the slope Sj in (4.6) implies that 
(4.7a) TV (R (. ; w)) = TV (w) ; 


hence it follows from (4.4) that the scheme is TVD, i.e., 


(4.7b) 


TV (E. • w) < TV (w) . 
h 



To simplify our presentation, we assume from now on that the partition 
is stationary and uniform, i.e. 


(4.8) 



5 


this enables us to express the schemes (4.3) by standard grid notations. 

The Godunov-type scheme (4.3) generates discrete values {v^^}, which are 
r-th order accurate approximations to cell-averages of the exact 
solution. We note, however, that the operation of the scheme (4.3) also 
involves a globally defined pointwise approximation to u(x,t) of the same 
order of accuracy which we denote by v^(x,t). The latter is defined for all 
X in the time-strips t^ < t < by 


(4.9) V ( •, t + t) - E(t) • R( • ; v^) for 0 < t <t 


We remark that (4.3) is the abstract operator expression of a scheme in 


the standard conservation form 
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(4.10a) 




with the numerical flux 


(4.10b) 




For r = 1 (Godunov's scheme), the numerical flux (4.10b) can be expressed 

by (1.7). For r > 2 , we make use of the fact that v, (x.. i, , t + t) in 

h j+ V2 

(4.10b) is needed only "in the small", in order to derive simple but adequate 

approximations to the numerical flux (see [16] for more details). 

We remark that (4.7a) is sufficient but not a necessary condition for the 

scheme E, to be TVD (4.7b). Other choices of the slope s, in (4.6), such 
n j 


(4.6b)' h • Sj = S (^j_l/W. 
or 

(4.6b)" h . Sj = m (2Aj_i/W , (Wj^.i-Wj_^)/2,2A^^ l^w) 

do not satisfy (4.7a); nevertheless the resulting scheme is TVD. This is due 
to the helping hand of the cell-averaging operator, wljich is not taken into 
account in (4.4), 

MUSCL-type schemes, as all other TVD schemes, are second-order accurate 
only in the Lj^-sense . In order to achieve higher-order of accuracy, we have 
to weaken our control over the possible increase in total variation due to the 
reconstruction step. We do so by introducing the notion of essentially non- 
oscillatory (ENO) schemes in the next section. 
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5. ENO SCHEMES. 

We turn now to describe the recently developed essentially non- 
oscillatory (ENO) schemes of [16], which can be made accurate to any finite 
order r . These are Godunov-type schemes (4.3) in which the reconstruction 
R(x ; w) , in addition to relations (4.2), also satisfies 

(5.1) TV (R( . ; w)) < TV (w) + 0(h^'’’P) , p > 0 

for any plecewlse-smooth function w(x). Such a reconstruction is essentially 
nonoscillatory in the sense that it may not have a Gibbs-like phenomenon at 
jump-discontinuities of w(x) , which involves the generation of 0(1) spurious 
oscillations (that are proportional to the size of the jump); it can, however, 
have small spurious oscillations which are produced in the smooth(er) part of 
w(x), and are usually of the size 0(h ) of the reconstruction error (4.2a). 

When we use an essentially non-oscillatory reconstruction in a Godunov- 
type scheme, it follows from (4.4) and (5.1) that the resulting scheme (4.3) 
is likewise essentially nonoscillatory (ENO) in the sense that for all 
piecewise-smooth functions w(x) 

(5.2) TV (E. • w) < TV (w) + 0(h^'''P) , p>0 ; 

n 

i.e., it is "almost TVD". Property (5.2) makes it reasonable to believe that 
at time t ® T, after applying the scheme N = T/t = 0(h ^) times, we can 
expect 

TV (v. (• , T)) < C • TV (u ) + 0(hP) . 
n o 


(5.3) 
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We recall that by Theorem 3.1, this implies that the scheme is convergent (at 
least in the sense of having convergent subsequences). This hope is supported 
by a very large number of numerical experiments. In order to conclude from 
(5.2) that for all n > 0 , 


(5.3) 


TV (v"'*’b < TV (v’^) + 0(h^‘‘‘P), p>0 


we still have to show that, starting from a piecewise-smooth 

u^(x) in (4.3b), v’^ remains sufficiently close in its regularity to a 
piecewise-smooth function, so that (5.2) applies to the following time-steps 
as well. Unfortunately, we have not been able as yet to analyze the 
regulatirty of v^ . 


Next we describe one of the techniques to obtain an ENO reconstruction. 



- ^4+ V2 

(5.4a) J w(y)dy = W(x^^ 1 ^^) - W(x^_ i^^) 

Va 

where 

(5*4b) W(x) - /^w(y)dy 

X 

o 


is the primitive function of w(x). Hence we can easily compute the point 
values {W(x^^ by summation 


(5.4c) 


W(x^^ 1/^ = h 


i 

I 


w. 


o 
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Let interpolation of u at the points 

accurate to order m, i.e. 

(5.5a) (yj ; u) = u(yj), 

Z i 

(5.5b) H (x * u) = — ^ u(x) + 0(h™^^ ^), 0<£.<m. 

dx dx 

We obtain our "reconstruction via primitive function" technique by 
defining 

(5.6) R (x; w) = -^ H^(x ; W). 

Relation (4.2a) follows immediately from (5.5b) with i = 1 and the 
definition (5.4) , i.e., 


R(x ; w) = -^ H^(x 5 W) = W(x) + 0(h’^) 


= w(x) + 0(h'^) 


Relation (4.2b) is a direct consequence of (5.5a) and (5.4), i.e.. 


V2 

A(Ij) R(. ; w) = i / ^ H^(x, W) 

"i- V2 


dx 


K = ”> - ■ K V2' - "<’‘1-1/2” 


W. 

J 
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To obtain an ENO reconstruction, we take in (5.6) to be the new ENO 
interpolation technique of the author 19], In this case, ; u) is a 
piecewise-polynomial function of x of degree m, which is defined (omitting the 
u dependence) by 


(5.7a) 


H_ (x ; u) = Y 4 <y<y 


m 


j ^ ^j+i 


where q. i. Is the unique polynomial of degree m that Interpolates u at the 

J+ V2 

m+1 points 


(5.7b) S^(l) 5 .... y^^) 

for a particular choice of i * i(j) (to be described in the following). To 
satisfy (5.5a), we need 


therefore, we limit our choice of i(j) to 


(5.7c) j - nri-1 < i(j) < j . 

The ENO Interpolation technique is nonlinear: At each interval 

[y^ , we consider the m possible choices of stencils (5.7b) subject to 

the restriction (5.7c), and assign to this interval the stencil in which u is 
"smoothest" in some sense; this is done by specifying i(j) in (5. ,7b). 
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The information about the smoothness of u can be extracted from a table 
of divided differences. The k-th divided difference of u 

(5.8a) , ...» = u[Sj^(i)] 

is defined inductively by 

(5.8b) u[S^(i)] » u(y^) 

and 

(5.8c) u[Sj^(i)J = (u(Sj^_j(i+l)] " u[Sj^_j(i)])/(yj^^j^-y^). 

If u(x) is m times differentiable in [y^^, ^i+m^ 

then 

(5.9a) u[S^(i)] “^“^”^(5) > for some y^ < 5 < 

If u^f*\x) has a jump discontinuity in [y^, ^i+m^ then 

(5.9b) u[S (i)] = 0(h““‘*’P[u^P^] ), 0 < p < m-1 

m 

( in the RHS of (5.9b) denotes the jump in the p-th derivative). 

Relations (5.9) show that |u[S^(i)]| is a measure of the smoothness of u 

in S (i), and therefore can serve as a tool to compare the relative 
m 

smoothness of u in various stencils. The simplest algorithm to assign 
S^(i(j)) to the interval [y^ . ^j+i- following: 
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Algorithm I . Choose i(j) so that 

(5.10) |„(S^ (1()))]1 - (|u 1S„(1)1|) . 

Clearly (5.10) selects the "smoothest" stencil, provided that h is 
sufficiently small (but not smaller than the round-off error of the machine 
would permit ! ) . 

In order to make a sensible selection of stencil also in the "pre- 
asymptotic" case, we prefer to use the following hierarchial algorithm: 


Algorithm II i Let be such that (ij^(j)) is our choice of a (k+1)- 

point stencil for Obviously we have to set 


(5.11a) 


(j) = j 


To choose consider as candidates the two stencils 


(5.11b) 


\+l * \+l ’ 


(5.11c) 


\+l \+l^^k^J^^ ’ 


which are obtained by adding a point to the left of (or to the right of) 
Sk^(ik(j)) , respectively# We select the one in which u is relatively 

smoother, i#e#. 


1^ (j) - 1 if lulS^ijll < |u[S^^jll 


(5. lid) 






otherwise# 



-31- 


Finally we set l(j) = 

Using Newton's form of Interpolation, we see that the polynomials 
{qj^(x)} , 1 < k < m , corresponding to the stencils selected 

by Algorithm II, satisfy the relation 

(5. lie) ^k+1^^^ “ FTfeCx-y) . 

yeS 

This shows that the choice made in (5. lid) selects q to be the one that 

deviates the least from q^^. It is this property that makes Algorithm II 
meaningful also for h in the pre-asymptotic range. 

In Figure 3, we apply the piecewise polynomial interpolation (5.7) to a 
plecewise-smooth function u which has in [-1,1] 3 jump discontinuities in the 
function itself and another one in the first derivative. This function is 
shown in Figure 3 by a continuous line on which there are 30 circles that 
denote the values used for the Interpolation. This function was continued 
periodically outside [-1,1] (not shown in the picture). 

In Figure 3a, we show the 6-th order polynomial (5.7) (i.e., m = 6) with 
the pfL2.deX2AmLn2.cL stencil l(j) = j; i.e., the 7-points stencil 

{yj. yj+i» •••. yj+ 5 > = S^(J) is used to define q^+ 1/^ ly^ . yj+il* 

Figure 3a shows a highly oscillatory behavior of the interpolation polynomial. 

In Figure 3b, we show the same 6-th order polynomial (5.7) except that 
now we use the adaptive, stencil which is selected by Algorithm II (5.11). 

To understand why this interpolation works as well as it does, we 
consider the following two possibilities: 

( i) [y^, is in the smooth part of u: For h sufficiently small, 

both Algorithms I and II choose a stencil S (i(j)) which is also in the 


smooth part of 


the function 


In this case. 


(5.5b) in 
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[y^, is the standard result for m-th order interpolation of a smooth 


function. We observe the q 


j+ V2 


need not be a monotone approximation to u in 


[y^, ^j+1^’ nevertheless, its total variation there cannot be more than 
OCh^+l) larger than that of u. 


(ii) [yj , Yj+j] contains a discontinuity: For h sufficiently small, 

the function u near [y^, can be thought of as a step-function. In the 

case of a step-function, the particular choice of i(j) is of no importance 

since all the stencils S (i) with j-m+1 < i < i lead to a q.. i, (x) which is 

m 3+ /2 

monotone in ^Yj » Yj^.j^)* This follows from the simple observation that in the 
case of a step-function, we have for all 1 < A < m, except A = j-i 


(5.12a) 




and, consequently, also 


(5.12b) 


V2^^i+ii ^ V 2 ^^i+A+1 ^ 


Using Rolle's theorem, we count in (5.12b) (m-1) roots of d q^^ j^^dx outside 
(yj , Yj^j) • Since dq^^ 1^^/dx is a polynomial of degree (m-1), it follows 
that these are all its roots. Hence, d q^^ j^^^dx does not vanish in 
(yj , Yj^j) » which shows that it is monotome there (see [17] and [15] for 
more details). 

We conclude this section by showing in Figures 4 and 5 the solution to 


the shock-tube problem (2.7) by the ENO scheme with r * 2 (Figure 4) and r = 4 
(Figure 5). Comparing Figures 4-5 to Figures 1-2, we observe a considerable 
improvement in performance (see [14] for more details). 
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6. NONLINEARITY, DPHIND DIFFERENCING AMD LINEAR STABILITY. 

In this section, we consider the constant coefficient case (2.1). In 
this case, the Godunov-type scheme (4.3) can be expressed as 

(6.1a) ’ 

where R(x; •) denotes the sliding average of R, i.e., 

, h/2 

(6.1b) R(x; •) / R(x+ 5 ; • ) d5 . 

-h/2 

We note that since R is a piecewise polynomial of degree (r-1), R is a 
piecewise-polynomial of degree r. Moreover, the conservation property (4.2b) 
shows that R(x ; v^) is an interpolation of interesting to 

note that using R which is obtained via the primitive function (5.6), we get 
from (6.1) the particularly simple form 


(6.2a) 


n+1 


= r [H (x 
h r 


j+ V 2 


- at; 


v”) - 


H (x. 1 , - ax; 


j- V 2 


V”) ], 


where defined at by (5.4), i.e.. 


(6.2b) 


n 


V2 


h E V. 

i»i ^ 


Relation(6 .2) directly relates Godunov-type schemes to interpolation. 
Clearly, if the interpolation is based on a fixed stencil, then the 
resulting scheme is linear; the nonlinearity of the ENO schemes stems from its 
adaptive selection of stencil. 
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When r = 1, R in (6.1b) is necessarily the piecewise-linear 
interpolation of » consequently the "upwind shift" (-ax) forces the 

scheme to be the first-order upwind scheme. We recall however that the 
stencil in the ENO scheme is chosen from considerations of smoothness which 
have nothing to do with the PDE; the "upwind shift" (-ax) is only by one 
cell; consequently the resulting ENO scheme (6.2) for r > 2 need not be, and 
in general is not, "upwind". 

We turn now to study the second order accurate ENO scheme (r = 2 in (6.1) 
- (6.2)). It is easy to see that this scheme is identical to the MUSCL-type 
scheme (4.6) with s^ defined by (4.6b)'. It is somewhat more surprising to 
find that the MUSCL-type scheme (in the constant coefficient case), is 
identical to the second order accurate modif ied-f lux scheme (3.11) with the 
correspondence 

(6.3) ^ * 

Consequently, all these second-order accurate TVD schemes can be written 
as (6.1) with a piecewise-parabolic R(x; v*^) . For a < 0, we get 

(6.4a) = R(xj - ax; v^) = Vj^ - X a jy-v*^ +V 2 a(l+Xa) • (s^^^-s^) , 

which is obtained from taking the sliding-average of (4.6a). 

We observe that when in (4.6b)' 



(6.4b) 
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then (6.4a) is the second-order accurate upwind-differencing scheme. However, 
when in (4.6b)^^ 


(6.4c) 


K ^ 

■ “i- Vj" ’ 


j+1 ^j+ ^2 


then (6.4a) is the central-differencing Lax-Wendroff scheme. 

Based on this observation, we see that the MUSCL-type scheme with 
defined by (4.6b) or (4.6b)^ satisfies (6.4b) when, as a function of i. 


(6.4b) 


{ 



is decreasing. 


and it satisfies (6.4c) when 


(6.4c) 



} is increasing. 


This shows that the ^popular” reference to the MUSCL scheme and the 
modified flux scheme as "upwind differencing" schemes is not justified. 

We remark that the scheme (6.4a) is second-order accurate only if 


(6.5a) 

This shows that in addition to 


“j+l ■ "xx • 


(6.5b) 


s. = u (x. ) + 0(h) 
J X J 


we need also the Lipschitz-continuity of the 0(h) term in (6.5b). As we have 
mentioned earlier, the MUSCL scheme, as well as the modified flux scheme and 
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Other TVD schemes, fail to have this extra smoothness at local extrema, which 
are the transition points between (6.4b) and (6.4c); consequently, their 
accuracy drops to first order at points of local extremum. 

The analysis of these second-order accurate nonlinear schemes shows that 
the "nature" of the scheme depends on differences of its numerical solution; 
therefore, local linearization is not justified. Since the two schemes in 
(6.4) are linearly stable, such incorrect linearization would nevertheless 
result in a correct statement of stability. This is not the case for 
r > 2, where, as r increases, more and more of the various choices of stencil 
can be identified as if belonging to a linearly unstable scheme. Since 

Fourier analysis is valid only if the same stencil is used everywhere, this 
identification is not necessarily relevant and may actually be quite 
misleading. 

A situation of this type is encountered when we consider the initial- 
boundary-value problem (IBVP) in -Kx<l 

(6.6) u. + u =0 , u(x, 0) = u (x) , u(-l,t) = g(t); 

t X o 

X = 1 is an "outflow boundary" and no condition needs to be specified there. 

We divide [-1,1] into (J+1) interval {I.}*? „ , where I. = (x. 

J 3=0 3 j 

and 



(6.7a) 



V 1/2 “ 


1 . 


Given cell averages 3 =0, ..., J we define W(x_ = 0 and 

compute j = •••> J by (5.4c) with i^ = 0; thus 1/^ I"® 
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given also at X = ±1. Next we evaluate H^(x; w) by Algorithm II, which is 
modified so that the choice of stencil in (5.11) is restricted to available 
data. Thus, H^(x;w) is defined for -Kx<l, and as before we define 

(6.6b) R (x; w) “ (x; w) , -Kx<l . 

Using this definition of R(x; w) in [-1,1] , we modify the Godunov-type 
scheme (4.3) by 

(6.7a) - A(Ij) • 'g(t) • R(« ; v") , 0 < j< J . 

(6.7b) Vj® “ * % * 

Here v (t) = E(t) • R( • ; v”) is the solution in the small (i.e., for 
0 < t < t) of the IBVP 

(6.7c) Vj. + “ 0, v(x,0) = R(x;v”) , v(-l,t) = g(t^ + t). 

This implementation of Godunov-type schemes to IBVP"s is very convenient: 
There are no "artificial numerical boundaries", and the prescribed boundary 
conditions are handled on the level of the PDE (6.7c). We observe, however, 
that near x = -1 the scheme is "differenced against the wind", which is 
linearly unstable if done everywhere. Therefore, our experience with linear 
schemes may inhibit us from using this approach. Overcoming this inhibition, 
we have performed a large numbe of numerical experiments with the so 
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modified ENO schemes (two of which are presented in the following) and we are 
happy to report that these schemes have been found to be stable in all our 
experiments. 

In Table 1, we present a mesh-refinement chart for the IBVP (6.6) with 
(6.8) u(x,0) = sin TTX , u(-l,t) = -sin tt (1+t) . 

The ENO schemes were used with a CFL number of 0.8, and the results are shown 
at t = 2. Table 1 indicates that the ENO schemes with 1 < r < 6 are 
convergent in this case; the accumulation of error seems to be linear. 
Comparing Table 1 to the periodic case (see [g]), we observe that the results 
for the IBVP are slightly better in the asymptotic range, which is to be 
expected. 

Next we consider the IBVP (6.6) with 

(6.9a) u(x,o) = e * , u(-l,t) = e^^*" , 

the solution to which is 

(6.9b) u(x,t) = e 

We observe that: (i) |u^^\x,t)| is a monotone decreasing function of x 

for all k and t. Consequently, if we apply Algorithm II to u(» ,t) we get 
i(j) = j in (5.11). (ii) The scheme (6.2) with the fixed choice i(j) = j is 
linear and strongly "biased against the wind"; consequently, it is linearly 


unstable 
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In Table 2, we present a mesh-ref ioement chart for the solution at t = 1 
of the IBVP (6.9) by the 4-th order ENO scheme (r = 4 in (6.6) - (6.7)) with 
CFL = 0.4. In spite of the previous observations, we find that the scheme 
seems to be convergent. This "paradox" is resolved once we examine the data 
in Figures 6 and 7 for J = 80 and 160, respectively. In (a), (b), (c) and 
(d), we show the absolute value of the k-th divided difference |v(Sj^]| for 
k = 0, 1, 2, 3, respectively. We see that the numerical solution and its 

first divided difference are monotone. However, the second and third divided 
differences are oscillatory. This allows the scheme to select i(j) ^ j in 
(5.11). The actual choice of i(j) at t = 1 is shown in Figure 6e and Figure 
7e; the straight line in these figures is i(j) = j. Comparing Figure 6d to 
Figure 7d, we see that the oscillations in vtS^] , the third divided 
difference of v, are uniformly bounded under refinement. Analysis of the 
numerical data suggests that 

(i) V [S^] — u^^\x,t) (in an average sense) as h -»• 0; 

(ii) v[S^] = u^’*^^(x,t) + for k=0, 1, 2. 

Finally, we consider the application of the 4— th order ENO scheme to the 
periodic IVP 

(6.10a) u + u =0, u(x,0) = ??? , v ° = (-1)^. 

L X J 


We observe that the mesh oscillation data in (6«i0a), 
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(6.10b) = (-l)J = 

is the highest frequency in (2.3b), which determines the linear stability of 
the constant coefficient scheme (2.2). We note, however, that as h decreases, 
the total variation of v° becomes unbounded. Consequently, v° does not 
represent a BV function and, therefore, need not be considered when testing 
for total-variation-stability in (3.1). In the following, we describe 
numerical experiments where we apply the 4-th order ENO scheme to (6.10) 
anyhow. The selection of stencil (5.11) is designed to make a sensible choice 
only when applied to piecewise-smooth data. In the mesh-oscillation case 
|v[Sj^(i)]| is constant as a function of i for all k; consequently, (5.11) 
results in the arbitrary choice of the uniform stencil i(j) “ j-3 (see Figure 
8b). As in the previous case, the ENO scheme becomes a constant coefficient 
scheme (2.2) for which linear stability analysis applies. In Figure 8c, we 
show the amplification factor of the mesh-oscillation mode 

k . 

(6.10c) g(v) = I (-1)*' C (v) 

i=-k 

as a function of the CFL number v = Xa. The amplification factor (6.10c) 
for the ENO schemes is determined by two competing factors: (i) Increase of 
oscillations due to the reconstruction, which is based on the highly- 
oscillatory interpolation of the mesh oscillation (6.10b); (ii) Decrease of 
oscillations due to the operation of cell-averaging on the translated data. 
Figure 8c shows that for v < 0.26, the latter wins and the scheme is linearly 
stable; for larger values of v the scheme is linearly unstable. In Figures 
8d and 8e, we show the numerical solution of the 4-th order ENO scheme with 
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V = 0.6 after a single time-step (n=l) and twenty time-steps (n=20), 
respectively. Clearly, the numerical solution blows up like (1.67)? 

It is amusing to realize that this "linear instability" is Itself 
"nonlinearly unstable" in the sense that any nonuniform perturbation of the 
mesh-oscillation data turns the ENO scheme into a stable nonlinear scheme. To 
demonstrate this point, we perturb the mesh-oscillation data by a random noise 
of the size 10 ^ of the round-off error (see Figures 9a and 9b), and repeat 
the previous calculation. In Figures 9d - 9k, we present subsequent "snap- 
shots" of the numerical solution, which show that the numerical solution 
decays in both the amplitude and the number of oscillations; observe that the 
rate of decay is faster for the highly oscillatory components of the solution 
and slower for the smoother ones. 

This property enables the scheme to combine "robustness" with accuracy. 
We demonstrate this feature of the ENO schemes in Figure 10 where we apply the 
4-th order scheme with v = 0.4 to initial data of sin irx perturbed by 
random noise of the size 10 ^ ; the squares denote the numerical solution; 
the continuous line shows sin irx. 
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Figure 5. Fourth order ENO scheme 
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Table 2 . Mesh RefinemenC for 4"th Order ENO with Exponential Data. 
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